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Abstract 


A technique is shown whereby it is possible to ielate a particular imiltignd pro- 
cess to cyclic reduction using purely mathematical aiguments This technique suggests 
methods for solving Poisson’s equation m 1-, 2-, or 3-dimensions with Diriclilet oi Neu- 
mann boundary conditions In one dimension the method is exact and, in fact reduces 
to cyclic reduction Tins provides a valuable lefeience point for understanding multigrid 
technique^ The particular multignd piocess analyzed is referied to here as Approximate 
Cyclic Reduction (ACR) and is one of a class known as Multignd Reduction methods m 
the literatuie It involves one approximation with a known enor term It is possible to 
relate the error term in this approximation with certain eigenvector components of the 
enoi These aie shaiply reduced in amplitude by classical relaxation techniques The 
approximation can thus be made a veiy good one 

1 lntioduct ion 


In the hist decade a new class of iclaxation schemes known as multignd methods have appealed 
in the literatuie Thes< scheme- solve laige, sparse, wide-banded linear equations They have 
nianv potential benefits including lemarkable convergence speed which i« -uch that the liiimbei of 
opeiations requited foi a solution is proportional to the liumbei of unknowns (ief 1 ) 

Multignd technique- appeal to be ideal for the solution of mail} equations emounteied in 
computational fluid dynamics today In addition to then speed they havt a potential foi simpler 
implement ation of component and solution adaptive grid- a- well as a liumbei of othei advantages 
detailed m (ref l) and (ref 2) Although widely expeiimented with, they do not appear to be 
in widespread use for practical pioblems The notable exception is the woik of Antony Jameson 
(ief 3) The purpose of this papei is to provide a clear and coheient explanation of a particulai 
multignd strategy when applied to a model pioblem and to show that multignd can be thought of 
as an approximation to cyclic reduction 

Multignd strategies ran most easily be analyzed when they are used to solve Poisson’s equa- 
tion This is because the matrix which approximates the Laplacian opcratoi has analytically known 
eigenvector and eigenvalue- Since the eigenvector aie sine functions, a Von Neumann stability 
analysis will usually agiee quite well with an analysis which use- (lie eigenvectors of the Laplacian as 
erroi components One exception occurs in the neighborhood of a boundary where the periodicity 
assumptions hi a Von Neumann analysis break down 

It is possible in one dimension to show a link between cyclic ieduction and multignd Although 
the analogy does not carry ovei exactly m two dimensions, it does provide some insight The 
analysis extends to two and three dimensions, and allows at least Neumann and Diriclilet boundary 
conditions 



Von Neumann stability analysis will usually agree quite well with an analysis which uses the 
eigenvectors of the Laplacian as error components One exception occurs in the neighbor- 
hood of a boundary where the periodicity assumptions in a Von Neumann analysis break 
down 

It is possible in one dimension to show a link between cyclic reduction and multigrid 
Although the analogy does not carry over exactly in two dimensions, il does provide some 
insight The analysis extends to two and three dimensions, and allows at least Neumann 
and Dirichlet boundary conditions 

Results are presented comparing the convergence rate of ACR with that of other clas- 
sical methods While the convergence of ACR is quite good compared to some methods its 
principal value lies in the insights it provides 

2 General Idea 


We assume some suitable discretization of the locally linearized governing equations 
Tins results in the system of linear equations 


~ ff 


(1) 


where Af is a matrix, <f>f is a vector of unknowns, and ff is a vector containing boundary 
conditions and a forcing function 

With multigrid techniques as with cyclic reduction the idea is to deduce the solution 
to equation (1) from the solution of a simpler equation 


Ac4*c — fr 


( 2 ) 


Cyclic reduction is able to exactly solve for 4>f from 6 C but mult.igrid can do this only 
in an approximate sense Tradit lonallv, each element of <f> c is approximately or exactly equal 
to a particular element of 4> j Also, for multigrid methods, A r is usually chosen to be of 
the same form as A j in sortie sense, which allows the solution of equation (2) to be derived 
from that of a still simpler equation In this paper we choose to examine only those schemes 
where A c and d> c are defined in this way 

Since both A c and <f> t have been chosen there must be a unique value of f c for which 
equation (2) is satisfied This paper will address the problem of how to find f c from the 
given information, A /, A c , and / / 

We begin by reordering the scalar equations and unknowns in equation (I) and per- 
forming the appropriate row and column permutations on A j 
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( 3 ) 


A\ 

1 A 2 

X 

4>t 


'fe' 

CO 

1 ^4 . 


<Po 


A 


This can be written as a set of two matrix equations with two vector unknowns 


A\4> t -t Az<f> 0 — f c 
A3<f> e -t Ai<f>„ — f„b 

If we solve equation (5) for <f> 0 we get 

4>o = A\ l (f 0 - A s <f> e ) 

Substituting this into equation (4) gives 

(>4i - A-2A 4 J A i)(f> e - ft AiA 4 f„ 
Notice that this has the form of equation (2) where 


( 4 ) 


(6) 


( 7 ) 


A r — (Aj - A?A 4 A f) 

(j) r (p e 

f c “ f e ~ A-iA 4 f„ 


(8a) 

(8b) 

(8c) 


By construction <f> t is a subset of <})f In general f c is not f e The nature of A r will depend 
on the nature of Aj and on which unknowns we choose to call <j> t We would like A c and 
Af to be identical discretizations of the same PDE on different sized meshes We will show 
that this is possible for Poisson’s equation 

3 One-Dimensional Example 


We will now focus on the particular set of equations which arise from discretizing 
Poisson’s equation with the standard central differencing In one dimension this is, for 
Dinchlet boundary conditions 
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'-2 1 


<t> i 


7i 

1 -2 1 


<t>2 


Si 


X 


- 


1 -2 1 


4>n - i 


Sn- 1 

1 -2 


_<j>N 


Sn 


(9) 


Here we have absorbed the boundary condition information into f\ and / jy 

In this case an elegant choice is to let <f> c be the even numbered dependent variables 
When this is done equation (3) becomes 


-2 

1 1 



<f>2 


'h 

.2 

1 

i 


4>i 


Sa 

-2 


i i 


<f>N-Z 


Sn - 3 

-2 


i i 


4>n - i 


Sn - i 

1 

_2 


X 

4> i 

— 

fi 

1 1 

-2 



<f>3 


fs 

1 1 


-2 


4>b 


/s 



-2 





1 1 


-2 


<t>N- 2 


fN-2 

1 


-2 


4>n 


Sn 


(10) 


For this example Aj 1 is just - 7/2 where 7 is the identity matrix Using this and multiplying 
by 27, equation (7) becomes 


'-2 1 


4> 2 


/l + 2/ 2 -f /3 

1 -2 1 


<t> 4 


/3 + 2/ 4 + Sh 


X 


r 


1 -2 1 


<t>N- 3 


Sn- 4 + 2 /. N -3 + / at- 2 

1 -2 


4>n - i _ 


f N -2 + 2//V-] + 


which is the desired result Note that A c in equation (11) has the same tridiagonal structure 
that A / had in equation (9) and therefore can be reduced in the same way If TV is one 
less than a power of two the reduction process can be continued recursively until only one 
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equation remains Once 4> e is known <f> 0 can be found by direct application of equation (6) 
This process is well known and is one of a class known as cyclic reduction (ref 4) It also 
represents a multigrid process where the restriction, sometimes known as the fine to coarse 
interpolation, is just 


(/c), - (//) 2 ,-i + 2 (//) 2 , + iff) 


2> + l 


and the coarse to fine grid interpolation stencil is just 


(<£/) 2t = [4>c) t 

/) 2 i -J- 1 = 2 ^ 1 ~ (■/'/) 2i+l 


where the required values at 0 and ( N + 1 )/2 are 


( 12 ) 


(13a) 

(13b) 


(0c) o — (‘/'tJjAH 1 )/2 — ^ (13°) 

In this way the restriction and interpolation are accomplished using the original difference 
equations By using the analysis techniques of cyclic reduction we are able to find an exact 
interpolation and restriction in one dimension 

4 Two Dimensions (Restriction and Interpolation) 


In two dimensions the multigrid processes defined by Brandt and others depart from 
standard cyclic reduction for the case of Poisson’s equation on a rectangle This comes 
mostly from the choice of <fi e 

If we index the unknowns as , corresponding to their x and y locations on the 
computational mesh, we see that cyclic reduction chooses <£ e to be those </>,j for which i is 
even The matrix A 4 is then block diagonal with each block a tridiagonal The inversion of 
such a matrix is just a series of one-dimensional problems, which is what makes it possible 
to compute f c In this case the matrix A c does not have the same form as Af but may be 
factored into a series of one-dimensional problems This approach is severely limited by the 
requirement that A c factor exactly 

On the other hand conventional multigrid, guided by physical intuition and a desire to 
reduce the number of unknowns faster, defines <f> c as those </>,,, for which t and j are both 
even For this choice all that can be said about A 4 in general is that it is wide banded 
and without any convenient structure An example will be shown in the next section The 
matrix A c , computed using equation (8a), is much more difficult to solve than Af and does 
not have the same form Some sort of approximation seems in order 
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One such approximation involves the standard decomposition 


A 4 — D + (L + {/) (14) 

were D,L , and U are diagonal, lower triangular, and upper triangular matrixes respectively 
Using this definition (5), (6), and (7) become 


^3<Ae+ \D + {L + U)}<j> 0 = f 0 ( 15 ) 

4o = D- 1 [f 0 -Aa4> t -{L+V)<l> 0 \ (16) 

(At - A 2 D- 1 A 3 )4>' = f e - A 2 D~ l f a + A 2 D~ 1 (L + U)(f> 0 (17) 

respectively 

It is often possible to express A 2 D~ 1 (L + U)<f> 0 in terms of (f> e and /,, That is 

A 2 D-'{L + U)4> 0 = G(f> e 4 H f a + t (18) 

where G and H are matrices which are chosen to minimize the error term and simplify A c 
The approximation is necessary for it allows us to eliminate (f> 0 from equation (17) It, is 
made possible by the fact that f t , (f> e , and (p n are related by a differential equation as well 
as by difference equations Substituting equation (18) into equation (17) gives 

(A, - A 2 D- 1 A , - G)4> e = f e - ( A 2 D - 1 - H)J 0 + c (19) 

As a convenience we may left multiply equation (19) by an arbitrary diagonal matrix 
D In this case the two-dimensional equivalents of equations (8) are 


A c - D[A X - A 2 D~ i A 3 - G) 

(20a) 

<£c = </*« 

(20b) 

f c = D[f'-(A 2 D-' - H)f u + t] 

(20c) 


Derivation of the matrices G , H , and D as well as the size of the resulting error term, 
c, depends on the particular problem and boundary conditions The coarse to fine grid 
interpolation involves approximations to the differential equation and is problem dependent 
The nature of these approximations is best illustrated by an example Such an example, 
that of a two dimensional Poisson problem is given here 
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5 Two Dimensional Example 


As a specific two-dimensional example let us take 


<f> XX + 4>yy = f{x,y) ( 21a ) 

with boundary conditions 

<t> (0,y) - /i (y) ( 21b ) 

<P*(C,y) - / 2 (y) ( 21c ) 

(x.O) - / 3 (x) (21d) 

4> (i,£) -- j 4 (x) (21e) 


discretized on a 4 / 4 equally spaced Cartesian grid so that At - Ay - = Again we absorb 
the boundary data and a factor of Ax 2 into the right hand side Also we adopt the double 
subscript notation where for example, j is expressed as <f > 2 2 

Tins discretization is 
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We permute this, choosing 4> c to be those <t> XJ for which t and j are both even 



Consequently equation (17) becomes 



where the term AiD~ l [L -|- U)<f> 0 is given by 


A 2 D-'{L + U)<t> 0 = 


- 2 { 4 > 1 \ + <^31 + <^13 + ^33) 

( + <^31 + ^ 33 ) 

— ( + $ 13 + ^ 33 ) 

~ ( + 2^33) 


We must relate this term to <t> c and f e , using the fact that <j> e is not independent of <f> 0 
The two are related both through the difference equations and through the underlying 
differential equation (21) Use of the former would lead to an exact solution but would not 
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be much cheaper than solving equation (22) directly Using the latter allows us to represent 
each line of equation (25) as a local discretization of equation (21) that is different from that 
used in equation (22) While this is cheap enough it will involve an approximation due to 
different truncation errors between the two discretizations of equation (21) We will show 
that it is possible to make the approximation a very good one by means of an appropriate 
relaxation process This is the nature of multigrid methods 

Looking at the first line of equation (25) we have 


_ 2 (011 + 4 > 31 + 4>13 + <£33) 


(26) 


Since we must use the differential equations to make our approximation we must know 
the physical location of these points in the domain We find them to be the four diagonal 
neighbors of the point (f > 22 If we expand each of them in x and y derivatives of (f> about the 
point cf > 22 we find 


~ 2 (^ n ^ ^ 31 ^ “*■ ^ 33 ) — — 2022 - 


Ax 


Ax V>+ — {<f> 

xxxx + b^iaryv + 4’yyyy) 


22 

(27) 


The numerical Laplacian / 2 2 may be related to the physical Laplacian A x 2 V 2 (/> at the point. 
(x 2 ,y 2 ) using the Taylor series analysis used in equation (26) This gives 


/ 22 


22 


(28) 


Adding equation (28) to equation (27) gives 


~ 2 (0J» + 03J + 013 + <^ 33 ) - — 022 - / 22 (^xzy v) 22 

One can make a similar argument for each line of equation (25) (see figure 1) 

We see that G and H are now fully determined since no terms involving (f> 0 remain The 
matrix D is then determined from the constraint that A c represents the same discretization 
of the problem as Af but on a coarser mesh In this case it is just 4 1 

For this example equation (19), after multiplication by D, is given by 
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-4 1 

2 -4 


-4 1 

2 -4 



4>n 



<t>41 


S' 

4> 24 



. <^44 . 



f 12 d / 21 4 /i3 ■+ /32 
/41 4 2/,j + /43 
/l4 + 2/23 -f /s4 
2/34 + 2/43 


which is indeed what we would like Notice that a restriction operator has been suggested 
by the mathematics This was arrived at by using all of the available information about 
the differential equations and the boundary conditions It is applicable at the boundaries 
as well as in the center of the grid Also of interest is that the expression for the error term 
is the same at the boundaries as it is in the interior of the grid This turns out to be very 
helpful 
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We now wish to recover <j> 0 from the now known <p e This may be done using the 
difference equations and an approximation If, for example, we wish to find <f> n we may 
use the now known (f> 22 , the boundary conditions, and the interpolation formula given in 
equation (29) centered around the point (1,1) instead of the point (2,2) Similarly, all of the 
unknowns with two odd subscripts may be found by this formula The remaining unknowns 
may now be found directly from the initial difference equations This interpolation strategy, 
using the original difference equations as much as possible, is characteristic of MGR methods 
in general 


6 Smoothing 


Ideally the error term 2A T 4 <f> xxyy in equation (30) would be zero In practice it seldom 
is In general there is very little that can be said a prion about this term One escape from 
this dilemma is to use the correction formulation of equation (l) For this we add Ajcp^ to 
both sides of the negative of equation (1) 


A - <t> j) - A ftfj - ff (31) 

Afe t f =r t f (32) 

The vector <f >^ is the current guess for <f>j (The more conventional notation, <^>", has 
not. been used here since n is used elsewhere ) The quantity r/ is called the residual and 
may be formed explicitly from known quantities Since <pf — <p l f - e^, solving for will 
yield <pf Also since equation (32) has the form of equation (!) all of the analysis developed 
for equation (l) will apply to equation (32) We adopt the notation that <f> c is a subset of 
e l j and f c comes from r)-, 1 c from equation (30), 

(/c)ll — ( r /)l 2 “t ( r /)21 + [r l f) 23 + ( r /)32 (32a) 

This notational convenience frees us from having to refer to the “error of the error” as the 
grids become successively coarser 

It is possible to solve equation ( 1 ) by means of the nonstationary Point-Jacobi relax- 
ation scheme (ref 5) 


4> t+1 = 4, 1 - h l (f - A# 1 ) (33) 

where h is a scalar iteration parameter For classical Point-Jacobi h — 1/4 

In this case it can be shown (ref 5 ), that the exact solution to equation (32) as a 
function of space and iteration number may be written as the double sum 
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M N i 

e ‘= EEIfo H A mn )c mrj X ran (34) 

m = 1 n = l r=l 

where for this problem (Neumann boundary conditions) 


A 

(-A mn)ij Sin 


M N 


mn 4 "f* 2 COS 

i(m - |)7T j(n - |)7T 


M 


sin 


N 


(34a) 

(34b) 


where as before i and j are space coordinates and t is the iteration number The coefficients 
c mn are determined from the initial guess This is simply a decomposition of the error into 
the eigenvectors of Af 

For this example the domain is a square of side £ Using the identities £ = A/Ax = 
N Ay, t — iAx, and y - jA y equation equation (34b) may be written 


X mn - sin 


(m - \)nx (n - \)ny 


sin 


We define the attenuation factor a mn as 


t 

G mn ~ |J (i T h A mn ) 

T = 1 


(34c) 


(35) 


We may now differentiate equation (34) directly to evaluate the error term of equation 
(30) when the process is applied to the correction equation equation (32) These are the 
same for each line of equation (32), namely 


2Ar 4 e Ixyv - 2 


2 N M 


s) EE*-' 


n= 1 m =- 1 


rn - 


n ~ 2 




Xr 


(36) 


which we refer to collectively as the error in f c We define the error in an element of (f> c as 
the difference between that element and the corresponding element of (j>f We may see the 
relation between the errors in / c and those in 4> c by looking at one eigenvector at a time 
The error in f c due to Xu is 
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8 \M) \nJ 


(37) 


c u&i iA'ii 


We observe that the error in f c is exactly an eigenvector of Aj interpolated to the 
coarse mesh It is also an eigenvector of A c We may therefore divide it by the corresponding 
coarse-grid eigenvalue to find the error in <j> c due to An This eigenvalue is exactly 


(■Mil 


(1 — |)2jt (1 - l)2n 

-4 + 2 cos — h 2 cos 


N 


M 


(38) 


Expanding the cosine terms about zero with a Taylor series and ignoring higher order 
terms gives 


( A c ) 1 1 = -((tf/M) 2 + (n/N) 2 ) 


(39) 


For this example M - N Dividing equation (37) by equation (39) gives the error in 
4> c due to An as 


~h(%) auc " x " 

Thus the eigenvector An is transferred to the coarse grid with second order accuracy 
In a similar manner one can show that for n < (A/2)-], m < {M/2) - 1 the error in 
<f> c due to the eigenvector A' mn is 


2(n/M) 2 (n/N) 2 \n - (l/2)] 2 [m - (1/2 )} 2 0 mn CmnXmn 
— 4 -f 2cos|(2n - 1 ) tt / iV] + 2cos[(2m - 1)7 r/M] 


(41) 


Of this group the worst case is when n ~ (A^/2) - 1, m = (Ai/2) - 1 If M and N 
are large this error approaches 


~ . G mnCmn^ mn 

b4 


(42) 


7'he eigenvector X mn shows the shape of the error The coefficient c mn depends on 
the initial guess The factor — 7r 4 /64 is roughly —15 Therefore, the attenuation factor a mn 
had better be less than 64/7T 4 in absolute value for this mode if it is to damp This is easily 
done As we will see in the next section, a much smaller value of a is required for certain 
other modes because of aliasing 
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7 Aliasing 


While there are MN fine grid eigenvectors there are only MN/4 coarse grid eigen- 
vectors Therefore some of the fine grid eigenvectors are not accurately represented on the 
coarse grid When m > M/2, or n > N / 2 or both we find that there is no eigenvector on 
the coarse grid corresponding to X mn on the fine grid The error e xxyy is not free from 
these components however, and when brought to the coarse grid, they appear as linear 
combinations of the coarse-grid eigenvectors More specifically, the error appears as the 
eigenvector X m ‘ n > on the coarse grid where 


m 


for m < M/2 


|M+l-m form>M/2, n 


n for n < N/2 

N + 1 — n for n > N/2 


(43) 


We now attempt to tailor a mn in such a way as to make a mn small for all values of m 
and n where X mn makes a large contribution to errors in <f> c Recall that 


A m „ - -4 T 2 cos 


( w ~ \)* , n ( n ~ \)* 


M 


+ 2 cos 


N 


(34a) 


and 


) ( 35 ) 

T — 1 

The parameters that determine o mn are h T For the moment, let us allow A m „ to 
have any value allowed by the range of the cosine terms This gives a two-dimensional 
space of A m „ which may be plotted as a square of side 2 centered at the origin The actual 
boundaries of the square are not included in the region (See figure 2 ) 

We see that on this diagram lines of constant A mn have a slope of -1 The attenuation 
factor (r mn will be constant along such lines as is evident from equation (35) Furthermore 
o mn will be zero when h T ~ ~1/X mn Any value of h < ~ A is strongly stable in the sense that 
kmn| < 1 for all m,n even if all the h T have this value One can also represent aliasing on 
this diagram by drawing contours of the fine grid A mn associated with eigenvectors which 
alias into coarse grid eigenvectors that all have t fie same value of A mn These contours are 
diamond shaped as shown in figure 3 

Before making a parameter choice, we briefly review the two major sources of er- 
ror First there is the term 2Ai 4 e Ixyv Equation (36) shows that this term is largest for 
the eigenvector X 'mm In this case, the term 2Ax‘ , e xiyy for small Ax, would evaluate to, 
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2n 4 a m n c m k X m n The shape of this mode is preserved but the amplitude is off by about 
a factor of 200 Actually this is only half of the story 

The second source of error is aliasing The largest aliasing errors occur on eigenvec- 
tors which alias into the coarse-gnd eigenvector associated with the smallest coarse-gnd 
eigenvalue From equation (39) this is, 


(A c )c-((7T/M) 2 + (n/N) 2 ) (44) 

Eigenvectors which alias into the eigenvector Xu on the coarse grid have their error 
term multiplied by the inverse of this quantity and deposited in X n The fine grid eigen- 
values in question are Xmn, Xmi, Xw, and Xu (The fine grid eigenvector Xu is in fact 
what we would like to transfer to the coarse grid eigenvector A'u It therefore can not be 
said to alias but is included here for comparison purposes ) Evaluating equation (41) by 
replacing m and n in the denominator with N' gives errors in <f> c of 


- 7r 2 N 2 O mn CmhXmn 

(44a) 

* 2 y 

T a 1 

4 

(44b) 

7T 2 

(44c) 

X\N 

4 

/ 7T \ 2 

(44d) 



respectively The first of these is unbounded if <tm/v is> 1 Clearly this can’t be tolerated 
We would like to pick h T such that om n 1s proportional to N' 4 By choosing M 1 ) — — | 

we are led to 0mn — (oiMjv) 4 Using equation (44a) this leads to an error on the coarse 
grid due to this mode of (1 /TV 2 )(7r G /64) which goes to zero in a second order way Thus we 
have completely neutralized the threat of aliasing from Xmn Errors from this source are 
of the same order as errors from Xu 

The errors represented by equation (44b) and equation (44c) must also be attenuated 
or they will dominate We are guided by a desire to preserve the accuracy of the eigenvector 
A’ji since this will be represented on even the coarsest mesh The relaxation selected in 
the previous paragraph also works on eigenvectors X\m and Ajvi, but not as well The 
attenuation factor for each step is = ctjvi = \ Since there were two steps, the total 
attenuation for these terms is one fourth The coefficients are then reduced from n 2 to 
7r 2 /4 We would like them to be proportional to N " 2 SO that errors in each of these modes 
decline in a second order way To do this we select h ^ ^ This results in a N} and <7j m 

both being equal to zero 

Finally, we look back to the errors m the components that don’t alias At the end 
of section 6 these were determined to be greatest in the eigenvector X 'm,n, which has a 
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corresponding eigenvalue approaching 4 This is the same eigenvalue as Xm i and Xm No 
further smoothing is required to reduce errors m components that don’t alias 

The question of smoothing parameters can be resolved for the model problem The 
coarse grid error component corresponding to the eigenvalue of smallest modulus is the one 
to protect All modes which alias onto this mode should be attenuated until the resulting 
coarse grid error decreases with N Their corresponding eigenvalues may often be estimated 
using Gerschgorin’s theorem (Ref 6) On problems which are within a perturbation of the 
model problem one might use the parameters given here scaled according to the largest 
eigenvalue 

For the three-step relaxation just discussed, (7 mn = [l + (A TOn /8)] 2 [l + (A mn /4)], which 
is plotted in figure 4 Thus, all of the eigenvectors which alias are severely reduced in 
amplitude and those which alias most are reduced most The eigenvectors which don’t alias 
are also reduced in amplitude 



8 Coarse To Fine Interpolation 


The errors incurred during the interpolation pose less of a problem than those incurred 
during the restriction process This is because the relevant error terms appear in the 
unknowns rather than in the right hand side Furthermore, aliasing does not occur since 
all the coarse grid eigenvectors are representable on the fine grid 

In what follows, we will explore and analyze one possible interpolation strategy This 
strategy is motivated by the same geometrical arguments used in forming the restriction 
operator 
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( 1 ) 


( e /)*J _ (^c)i i 


if * and j are even 


( 2 ) 


( e /)*j ~ ^ l( e /)t+Xj + l + ( e /)»-lj+l 

+ (e/) t +ij-i + (e/),-ij-i 

"2 (r/)„ 1 


if « and j are odd 

( e /)*x ~ ~ 4 I( e /)t+ij + i + ( e /)«j + 1 

+ ( e /)«~ij + ( e /)tj-i 

- ( r f)n ] 

for all t and j This comes directly form the original difference equations for these 
points 

(4) Improve the estimate of <f>j by subtacting ej 

An explanation of the above strategy follows In step 1 we simply assign coarse grid 
values to the fine grid at those points where the two grids coincide In step 2, we again use 
the rotated difference equations used during the derivation of the restriction operator This 
carries with it a fourth order error term which causes inaccuracies in all of the eigenvector 
components of < j In step 3, we use the difference equations to fill in all the missing values 
Finally, in step 4, we use our knowledge of the fine grid error to improve the fine grid 
solution Optionally one can do some more smoothing to remove the errors incurred during 
this interpolation process Though this improves the convergence per step it was found not 
to be cost effective See the discussion of operation count for more details 

9 A Three-Dimensional Example 


(3) 


The example chosen in three dimensions is again Poisson’s equation using purely 
Dirichlet boundary conditions Again we use an equally spaced Cartesian mesh and the 
standard seven-point differencing star The domain is a cube of side £ The analysis is 
exactly the same as for the two-dimensional problem although the matrices are much larger 
Only the results of the analysis will be given here 

We are led to the following approximation in arriving at the restriction operator 
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“($t+l 3 k+l+0i-l j Jfc +1 d~ <£. + 1 3 k-i d- 3 fc -1 

+ 0t 3 + I fc+l + ^i 3 + 1 fc-J + 0i 3 -I fc+1 + 0i 3 -I fc-1 

+ <f>i + l 3 + 1 k + 0 i-l 3 + 1 fc + <£.+ 1 j-1 fc + <^t-l j-1 fc) 
A 

— 301] fc+ /,j fc + ~ ( 0xxyy d <t>yvzz + 0 zzxx) 


This yields the restriction operator 

{fc)tjk= (17)2.41 23 2fc + ( r /) 2. — 1 23 2 k + (r 7)2. 23 + 1 2fc 

+ i r f)2% 23-I Ik + ( r /)2. 23 2fc + l + ( r /)2. 2j2k-l - 2 (r/) 2 , 23 2fc 

In analyzing the smoothing, the three-dimensional analog of figure 2 is a cube instead of a 
square Reasoning, as in section 7, we choose four relaxation sweeps with values of h T = 
j3j , |, and d The last is not strongly stable but the sequence is stable Notice the similarity 
to the two-dimensional case where the worst error is smoothed twice and the other errors 
are smoothed once 

The coarse to fine-grid interpolation is similar to that for two dimensions although one 
more approximation is required 

1) 

for i,j,k, even 


( 2 ) ( e /).3fc — g K e /)«4l34 1 fc+l d ( c /). + l34 lfc -1 

+ ( f /). 4 13 — 1 fc4 1 + ( e /)«4 lj-lk 1 

d ( e /)« + l 3 +lfc+l d- ( e /).4l34 lk - 1 
d ( e /).4 I 3 -U + I d (e/),+ 13 - lk 1 
-4{rf) ljk j 


if 1 , j, and k are odd 
(3) ( c f)ijk = 


^ l( e /)»+l3 fc d" { e f)i-ljk } 
d-g [( e /).3 + U41 d- (e/),j + lfc_l 
+ ( e /)«3 -U+l d- (e/)., - 1 fc— 1 ] 
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for t odd, j and k even and for k off, t, and j even 


( 4 ) 


( e f)tjk ~ g [ { e f)t + l]k + ( e /)x-ljfc 

+ ( e /)*J + lfc + i e f)i]~lk 

+ + (e/),jfc_l 

~ ( r /).jfc ] 


for all indices This is just using the original difference equations 
5) Improve the estimate of 4>f by subtracting ey 

10 Summary 


The process just described is summarized as follows 

1) Smoothing to reduce the error incurred during restricton This is done using nonstation- 
ary Point Jacobi relaxation with the h l selected above 

2) Computation of the required fine grid residuals 

3) Transferring the problem to a coarser mesh using the restriction operator derived above 

4) Exact solution of the problem on the coarser mesh If the coarsest mesh has more 
than 1 unknown, “exact” solution may be the result of some suitable relaxation process 
This will be cheap since the coarsest mesh has very few unknowns On other than the 
coarsest mesh “exact” solution means two iterations of this multigrid process (This is 
the so called W-cycle) 

5) Transferring the solution back to the fine mesh using the coarse to fine interpolation 
given above 

6) Repetition of steps 1 through 5 until convergence is obtained There will be further 
discussion of what, is meant by convergence 

11 Operation Count 


In this section we address the total cost of ACR In two dimensions, a nonstationary Point 
Jacobi relaxation for the five point Laplacian requires 7 operations per point where mul- 
tiplications and additions are both counted Interpolations account for about 30% of the 



Table 1 ACR Operation Couni 

2-D 

Description Of Process Segment 

3-D 

21 

finest mesh smoothing 

36 

3 

computation of residual? 

4 

1 

restriction operator 

7 

8 

4 

interpolation operator 


2') 

total for finest mesh 

47 1 

y 2 

factor for W-cycle 

x| 

58 

total for all meshes 

63 f 


total For three dimensions the relaxation sweeps require 9 operations/point Interpolations 
require only about 25% of the total The operation count for both is given in table 1 

In both of the above operation counts we have taken into account the fact that restriction 
only occurs at fine mesh points with all even subscripts This means that we do not need 
the residuals everywhere The count for both of these reflects the fact, that they are not 
done at every point The factor for the W-cycle assumes an infinite number of grids In the 
two-dimensional case for example, each grid requires one-fourth the number of operations 
of the next finer grid, but must be visited twice for each time the finer grid is visited This 
leads to the series 


1 + 


1 

2 



1 

+ 8 -* 


- 2 


which is where that factor of two comes from In three dimensions each grid requires only 
one-eighth the number of operations of the next finer grid This leads to a factor of four 
thirds Notice, that because the number of operations on coarse grids is proportionately 
less in three dimensions than in two dimensions, the cost of an additional relaxation sweep 
is also less even though the difference stencil is larger 

This scheme was devised for ease of explanation rather than for speed Possible speed 
improvements include 

1) Improved relaxation schemes such as checkerboard Gauss-Seidel or incomplete LU de- 
composition These schemes are more efficient at removing all the restriction errors and 
require no parameter choice, but are more difficult to analyze 

2) Configuring the scheme as an FMG cycle (ref 1) In this case, the method would start 
with an exact solution on the coarsest mesh It would then proceed as described above 
but starting at the point where the coarsest grid exact solution is computed This can be 
thought of as producing a better initial guess on the finest grid at minimal cost Some 
investigators have found that only one additional cycle is required to reduce the errors 
to the level allowed by our finite difference approximation It rarely makes sense to 
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reduce the errors to a still lower level All the cases run here were converged to machine 
accuracy however, since our intent was to test the convergence properties of the method 

3) Visiting each mesh only once instead of twice This is the so called V-cycle This gives a 
savings of one-third in the operation count per step (one-seventh m 3-D) but weakens the 
bounds on the spectral norm of the method because the ’exact’ solution on intermediate 
grids is not as good Consequently more steps may be required There is some practical 
experience to the contrary (ref 7) 

12 Invariant Subspace Analysis 


The question of errors introduced during the interpolation was only touched on briefly In 
fact it is these errors that allow reintroduction of high-frequency error on the finest grid 
Without them the troublesome components would soon disappear and the restriction would 
become nearly exact With exact interpolation and restriction multigrid becomes a direct 
method Since these errors limit convergence, it is necessary t,o take them into account when 
analyzing rnultigrid methods The best (perhaps only) quantitative analysis of interpolation 
errors for the model problem is the method of invariant subspaces This is explained in some 
detail in (ref 7) 

lri the section on restriction it, was sliowm that in two dimensions the four line grid eigenvec- 
tors A’ mn , A 'm'n, A'mn'i an d A” m ‘ all appear on the coarse grid as the A’ m „ eigenvector In 
the section on interpolation we briefly outlined a way of treating errors that only occur on 
certain points If we pursue this, we find that when the coarse grid eigenvector X mn is inter- 
polated to the fine grid, errors are introduced in only the four eigenvectors just mentioned 
This nice property is preserved through the smoothing restriction, and computation of 
residuals as well Thus the error in these four components at the end of a step depends 
only on their errors at the begimng of that step We can analytic ally form the 4 > 4 matrix 
which represents this situation To find the error at the end of T steps we simply multiply 
the initial error by the Tth power of this matrix The spec tral radius of the method is just 
the largest of the spectral radii of these 4x4 matrices and t ho sped ral norm is the largest 
of their spectral norms bsmg this method we have numerically computed these quantities 
for the methods advocated here Also we computed these quantities for different amounts 
of smoothing Adding a smoothing sweep can decrease the spectral radius and norm but 
will increase the cost A function w'hich follows this tradeoff is 

log of spectral norm 
number of operations 

Loosely, this is the number of base e digits per multiply 

In tables 2 and 3 we show the performance of ACR with different amounts of smoothing 
For completeness we include the possibility of post interpolation smoothing The numbers 
are limiting values as Ax — > 0 

The three-dimensional case is completely analogous although the invariant subspaces each 
contain eight components instead of four Consequently there are three eigenvalues cor- 
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Table2 2-D Convergence Results 

Smoothing 

Spectral 

Spectral 

F 

t 

Pre- Restrict ion 

Post-Interp 

Radius 

Norm 


0 






1 255 

oo 

— OO 

1 

i 

8 





0 500 

1 250 

-0 07 

2 

1 

8 

1 

8 




0 249 

0 390 

0 21 

3 

1 

8 

1 

8 

1 

4 



0 114 

0 118 

0 37 

4 

1 

8 

1 

8 

1 

4 

6 


0 074 

0 092 

0 33 

5 

* 

1 

8 

1 

4 

1 

6 

2 

U 072 

0 073 

0 30 


Table 3 3-D Convergence Result 5 ; 


Smoothing 

Spectral 

Radius 

Spectral 

Norm 

F 

t 

Pre-Restriction 

Post-Inlerp 

0 







3 323 

oo 

— OO 

1 

12 






1 188 

OO 

— oo 

2 

1 

12 

1 

12 





ii 297 

1 001 

-0 000 

3 

1 

12 

1 

12 

1 

8 




0 22U 

0 302 

0 017 

4 

1 

12 

1 

12 

1 

8 

\ 

4 



0 148 

0 192 

0 023 

5 

1 

12 

1 

12 

1 

8 

x 

4 

1 

c 


0 111. 

O 100 

0 022 

G 

1 

12 

1 

12 

x 

8 

1 

4 

1 

6 

1 

fl 089 

ll 132 

U U21 


responding to error components aliasing into the smallest eigenvalue instead of two (and 
seven eigenvectors instead of three) 

13 Results 


The 2-D results are for the example in section 5 Eight mesh sizes varying between N = 2 
and N — 256 were tried All the test cases were reducible to one unknown We chose the 
homogeneous case where f j <f> f — 0 This was chosen to simplify computation of the 
error (which for this case is just the current estimate for <f>f) and does not imply that the 
process is restricted to homogeneous boundary conditions (ref 8) 'The initial guess was 
chosen so that all the coefficients c mn were equal and of such a magnitude as to make the 
L 2 norm of the error equal to 1 The complete convergence history is given in table 4, for 
a 256 X 256 grid using the three step smoothing suggested in section 7 

Each complete cycle reduced the L 2 norm of the error by a factor of 27 or more The 
spectral norm guarantees a factor of 8 5 per step but this is overly pessimistic Like any 
linear iterative scheme, ACR starts out fast and then slows down to some convergence rate 
which depends on the spectral norm or radius Us advantage is twofold, first, the asymptotic 
rate is independent of N , quicker than any other explicit method and second, the problem 
may well be converged before this limit is reached The independence of the spectral norm 
on N has been shown for other multigrid methods (ref 7) Our experience has shown no 
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Table 4 2-D Convergence History — 256 x 256 Grid 

Step 

II « lb 

II e IL 

0 

] x 10° 

5 3 x 10 4 

1 

1 3 x 10' 2 

3 1 x 10 2 

2 

2 3 x 10‘ 4 

6 9 x 10° 

3 

2 6 x 10' 6 

5 0 x 1CT 2 

4 

6 3 x 10' 8 

7 9 x 10* 4 

5 

2 3 x lO' 9 

4 6 x 1CT 5 


degradation of convergence over a wide range of values for N 

The 3-D results are for the example in section 9 Four mesh sizes varying between N — 3 
and N — 31 were tried All the test cases were reducible to one unknown We again chose 
the homogeneous case where {/ - <f>f — 0 The initial guess was chosen so that all the 
coefficients c mn0 were equal and of such a magnitude as to make the Li norm of the error 
equal to 1 The complete 3-D convergence history is given in table 5 for a 31 x 31 x 31 grid 
using the four step relaxation suggested m section 9 


Table 5 3D Convergence Hi«tory-31 x 31 x 31 Grid 

Stej) 

II e lb 

II ' IL 

0 

1 0 x 10° 

2 3 X 10 4 

1 

2 ft x irr 2 

3 n x io 2 

2 

6 6 x 10' 4 

9 2 x 10" 

3 

2 6 x nr 8 

3 4 x 10 _! 

4 

2 9 x nr 0 

1 h x 10' 2 

5 

7 3 x 10' 8 

8 5 x 10‘ 4 


Each complete cycle reduced the Li norm of the error by a factor of 9 or more The spectral 
norm guarantees a factor of 5 2 per step 

Thus we see that the remarkable results claimed by the analysis are actually realized in 
practice No other type of explicit method allows an entire convergence history of this 
problem to be written down in a short table Furthermore, 3-D problems take only about 
1 5 times as much work per point as 2-D problems, an important feature of multigrid 
methods 

In comparing ACR with other explicit, methods and with cyclic reduction I will use table 
G prepared by Dorr (ref 9) This is for Poisson’s equation discretized on a square with 
N 2 unknowns The direct methods are compared with the iterative ones by assuming that 
a reduction of the error by a factor of N 2 is required 'This comes from the fact that 
the truncation error is proportional to jji For purposes of comparison all acceleration 
parameters are optimally chosen 

The factor of log 2 N in ACR and MGR-CHi.i does not come from any specific feature 
of the algorithm but from the fact that the desired accuracy increases with the number 
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Table 6 Method Comparison 

Method 

Operation Count 

Block (polynomial form) 

C,N ? 

Block (Schecter) 

| N' 

Block (Froehlich) 

(f‘ -+ 4r)N 3 

Odd-even reduction (Buzbee et al ) 

j log 2 N 

Tensor product (Lynch et al ) 

8N J 

Fourier Series (Hockney) 

f>A' 2 log 2 A 

SOR 

-j A' 3 log ^ N 

ADI 

4 A’ 2 (log , A') 2 

ACR 

~ 3s A' 2 log 2 N 

MGR-CHj,i 

~ I7A 2 log 2 N 


of unknowns Limited precision on a given computer may limit the attainable accuracy 
Under such conditions the method requires order iV 2 operations to achieve this limited 
precision In any event-, given the rest rictions on memory size common in today's comput ers, 
log 2 N 10 In practice multigrid methods can be made to be of order A' 2 by use of the 
FMG cycle outlined previously in reference 1 but the coefficients will increase from 38 and 
17 to 96 and 50 for ACR and MGR-ClU j respectively The FMG cvelo is usually good if 
the initial guess is largely random If, on the other hand there is a reasonable guess from 
some nearby problem the basic W -cycle will probably converge in one or two steps 

The improved performance of MGR-CH 2 J over ACR is due largely to its use of checkerboard 
Gauss-Seidel for the removal of high frequency error components The operation count for 
this method is much less than for the non-stationary point -Jacobi relaxation used in ACR 

14 Conclusions 


For the multigrid process just, presented it is possible to formally analv/e errors made in the 
interpolation and restriction processes on these model problems It then becomes possible 
to tailor the smoothing according to these errors The* analysis yields interpolations and 
restrictions that are valid at Neumann boundaries as well as in the interior of the domain 

We have shown that cyclic reduction can be thought of as a particular multigrid method that 
has exact interpolation and restriction This is particularly evident in one dimension where 
the two methods coincide Although the two methods differ in higher dimensions, they are 
equivalent up to a known approximation This approximation can be improved with an 
appropriate relaxation Using the difference equations to do the interpolation eliminates 
the need for post-interpolation smoothing An efficient explicit method results 

The real value of multigrid techniques comes from applications to problems which cannot 
be solved with cyclic reduction Since ACR can be viewed in terms of point operators it 
may prove easier to adapt to complicated grid structures than cyclic reduction 
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